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Abstract. We show that a cooled region of shocked supernova ejecta forms in a type II supernova-QSO wind 
interaction, and has a density, an ionization parameter, and a column density compatible with those inferred for 
the high ionization component of the broad emission line regions in QSOs. The calculations are based on the 



assumption that the ejecta flow is described initially by a similarity solution investigated by Chevalier ( 1982 ) and 



Nadyozhin ( 1981: ) and is spherically symmetric. Heating and cooling appropriate for gas irradiated by a nearby 
powerful continuum source is included in our model, together with reasonable assumptions for the properties of 
the QSO wind. The model results are also in agreement with observational correlations and imply reasonable 
supernova rates. 
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1. Introduction 

In active galactic nuclei (AGN), the primary physical 
mechanism for the excitation of the broad emission line 
regions (BELR) is photoionization by the underlying 
broad-band contin uum (e.g. Osterbrock & Matthews 1986| ; 
Clavel et al. 1991 ). Evidence for at least a two-component 
structure of th e BELR has been provided by e.g. C ollin- 
I1986Q and Wills et al. (p.985|). 



ii) cloud formation in shocks produced by the interaction 
of an accretion disc wind with a nuclear wind (Smith & 
Raine |1985| ); iii) the interaction of an outflowing wind 
with the surface of an accretion disc (Cassidy & Raine 
1996); iv) interaction of stars with accretion discs (Zurek 



et al. 1994); v) tidal disruption of stars in the gravitational 
field of the BH (Roos |1992fl; vi) interaction of an AGN 



Souffrin et al. ( |1982| , |1986j ) and Wills et al. (p.985[) . One 
component can be identified as consisting of high ioniza- 
tion lines, including Lycv, Cm], Civ, Hei, Hen, Nv, and 
other multiply ionized species, and is known as the HIL. 
The second component can be identified with the low ion- 
ization lines which include the bulk of the Balmcr lines, 
and lines of singly ionized species (e.g. Mgn, Cn and Fen), 
and is known as the LIL. The regions emitting the LIL and 
HIL display different kinematics, as deduced from studies 
of the profiles and line widths (e.g. Gaskell 1988, Sulentic 
et al. 1995| ) ■ ^ n quasars, the HIL are also systematically 
blue-shifted with respect to the LIL (see various papers in 
Gaskell et al. 199E). Discussion of the BELR parameters 



for a range of AGN properties can also be found in the 
articles in Gaskell et al. 



wind with supernovae an d star clusters (Perry & Dyson 
1985 ; Williams & Perry 1994 ); an d vii ) emission from 
accretion shocks (Fromerth & Melia 2001 ) . Models identi- 
fied as containing serious difficulties include the formation 
of BELRs by the thermal instability of a hot opti cally 
thin flow (as e.g. d iscussed by Beltrametti 1981 and 
Shlosman et al. |l985|) . Perry & Dyson ( |l985| ) noted that 
this mechanism will not occur when Lf, i > 10 46 ergs _1 
as Compton cooling dominates over bremsstrahlung, and 
therefore it cannot be responsible for the observed BELR 
in high luminosity QSOs. Two-phase equilibrium models 
(Krolik et al. 1981) are also not applicable to BELRs 
in these sources because implausibly high values of the 
AGN mass-loss rate would be required (Perry & Dyson 



1999 



Many theoretical explanations have been proposed for 
the origin of the BELR. They include: i) magnetic acceler- 
ation of clouds off accretion discs (Emmering et al. 1992|); 



1985). The formation of the BELR in ionized red giant or 
supergiant winds (Scoville & Norman |1 988 Kazanas 1989; 
Alexander & Netzer 1994) has difficulty in reproducing 
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the observed broad line wings (Alexander & Netzer 1997), 
and models involving the ballistic deceleration of clouds 
have a number of problems, as summarized by Osterbrock 
& Matthews (1986). There are also concerns about the 
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various assumptions in the infalling and orbiting cloud 
models proposed by Kwan and Carroll (Kwan & Carroll 



1982; Carroll 1985; Carroll & Kwan 1985) 



One model which accounts for the LIL emission was 
proposed by Collin-Souffrin et al. ( 1988Q . In this model 
the LIL emission arises from the surface of the accre- 
tion disc, which is illuminated by back-scattered X-rays 
from the central source. Typical electron densities and 
absorption columns are inferred to be n e > 10 11 cm -3 
and Nh Si 10 24 cm -2 respectively. Approximately three 
quarters of the total luminosity of the broad-line emis- 
sio n is e stimated to arise in this fashion (Collin-Souffrin et 
al. 1988J ). The HIL, therefore, contribute about one quarter 
of this emission. The characteristic mass of the B ELR i n 
high luminosity QSOs is - 100 M Q (Osterbrock |l993| ), 
and this is dominated by the HIL. 

In this paper we look at one aspect of activity in AGN 
connected with stars, and which is relevant to the work 



of Perry & Dyson ( 1985 ). We investigate the interaction 
of supernova ejecta with the optically thin, low density, 
hot QSO wind, in the presence of intense continuum radi- 
ation. In particular we examine if shocked gas can radiate 
efficiently enough to cool to temperatures appropriate for 
the HIL. The evolution of SNRs in a high density static 
ambient medium has been previously studied by Terlevich 
et al. (1992), with particular application to the formation 
of BELRs in starburst models developed to obviate the 
existence of supermassive black holes in AGNs. Although 
there are similarities between this work and ours, two 
differences exist. First, the initial conditions for the super- 
nova ejecta differ from those in our work. Second, these 
authors did not include Compton cooling or any heating 
processes in their calculations. These factors will influence 
the thermal evolution of the shocked regions. 

In Section 2 we discuss the use of a similarity solution 
to specify initial conditions in the calculation and the cool- 
ing and heating rates adopted. Section 3 contains results 
for the calculated evolution of a remnant for each of sev- 
eral assumed sets of environmental conditions, showing 
that the formation of a cool region of shocked ejecta hav- 
ing a density, ionization parameter, and column density in 
harmony with those inferred from observations occurs for 
reasonable assumptions. In Section 4 we summarize our 
conclusions and describe how the work can be extended. 

2. Details of the calculations 

Fits to the results of explosion models of type II su- 
pernovae indicate that power-law stratifications represent 
adequate approximations to ejecta density and velocity 
profiles (e.g. see Woosley et al. 1988 and earlier work 



by Chevalier 1976 and Jones et al. 1981), and have been 
widely used in analytical and numerical studies of remnant 
evolution. If p oc r~ n , these results indicate that n ss 12 
for the high velocity ejecta. 

Self-similar solutions for the structure of the shocked 
ejecta and swept-up medium assuming power-law approxi- 
mations for both of the unshocked components are derived 



by Chevalier (1982) and Nadyozhin (1985), following ear- 
lier work by Parker (1963). Their relevance to actual 
remnants was most recently highlighted by SN 1993 J 
in M81. High resolution, spatially resolved VLBI obser- 
vations showed self-similar evolution of the azimuthally 



averaged radius (Marcaide et al. 1997). From application 
of the Chevalier-Nadyozhin model, the observations are 
best fitted with n ~ 12, in good agreement with results 
from explosion models. 

For n > 3, there must be an inner core of material 
with a shallower density profile (p oc r ) in order for the 
mass of the ejecta to be finite. Such a core can be seen 
in the results of e xplosi on models (e.g. Jones et al. 1981 , 
Suzuki & Nomoto 1995 ). In the simplest case, which we 
adopt in this paper, a core with uniform density (5 = 0) is 
surrounded by a steep outer envelope (n = 12). The speed 
of the core radius i> c , and the density of the envelope p e (r) 
are then given by 



2(5- 6)(n-5)E 



(3-8){n-3)M ej 



1/2 



Pe = g n t~ 3 



(1) 

(2) 



where E is the explosion energy, M e j is the ejecta mass, 
and 



1 [2(5-J)(n-5)£]("- 3 )/ 2 

47r(n - S) [(3 - 6){n - 3)Af ej ](™- 5 )/ 2 



(3) 



(cf. Chevalier & Fransson [1994J). For 5 = and n = 12 
(which we assume in this paper), the density of the core 
p c is given by 



Pc 



729 



3 Ml 



1/2 



1120tt \ 70 E 3 



(4) 



(cf. Band & Laing 1988), whilst for <5 = and variable 
n the relative mass and energy of ejecta in the envelo pe 
compared to the core are (cf. Truelove & McKee 1999 ): 



M r 



M t 



Ez 

E 



n — 3 

n 
n — 5 

n 



(5) 
(6) 



Hence for n — 12, 75 per cent of the mass and 58 per cent 
of the explosion energy are in the core. 

If the explosion occurred in a pure vacuum, the struc- 
ture of the ejecta would evolve according to Eqs. [l]-|| 
However, when the ejecta interact with a surrounding 
medium, the steep envelope acts as a compressible piston, 
and the radius of the core relative to the reverse shock in- 
creases with time. As long as ejecta in the steep envelope 
continue to pass through the reverse shock the solution 
is self-similar, but once the core radius reaches this point 
the solution abruptly ceases this behaviour. 
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Fig. 1. The density profile for ejecta from a type II SN 
explosion, modelled as an n = 12 power-law for the en- 
velope and a 6 — power-law in the core, an explosion 
energy of 10 51 erg and an ejecta mass of 10 Mq. The dot- 
ted line shows the solution given by Eqs. [jlS which is valid 
for ejecta expanding into a total vacuum. The solid line 
shows the case when the progenitor is surrounded by a 
constant density medium with n — 10 4 cm -3 : the outer 
part of the ejecta and the swept-up ambient medium are 
compressed into the Chevalier-Nadyozhin similarity form. 



In this work we adopt the Chevalier-Nadyozhin 
similarity solutions to specify our initial conditions, with 
the ejecta distribution being specified by an inner core 
with (5 = and an outer envelope with n = 12. Whilst this 
is a simplification to actual distributions obtained from ex- 
plosion models, at this stage we aim to keep our results as 
general as possible. It also has the benefit of being easily 
reproducable and scale-free. In all our calculations we as- 
sume a canonical explosion energy of 10 51 ergs and ejecta 
mass of 10 Mq, which is typical of a SN of type II. Fig. [y 
shows a typical density profile at t = 0.1 yr. 

Heating and cooling rates for a canonical AGN spec- 
tr um w ere kindly supplied by Tod Woods (cf. Woods et 
al. 1996) and are included in our calculation. We made 
use of an adap tive grid hy drodynamic code (see e.g. Falle 
& Komissarov 1996) , 199§| ), which is ideally suited to this 
problem where regions which contain small-scale structure 
are located within much smoother regions. 

To test the accuracy of the code we first imposed the 
similarity structure on the flow and saw whether it sus- 
tained itself. The ambient medium had constant density, 
zero velocity, and negligible pressure. In Fig. |2| we show 
the results of this test, which compare favourably with 



the results from other codes in the recent literature (e.g. 
Blondin et al. |2000| ). 

The heating and cooling rates are tabulated as func- 
tions of temperature and ionization parameter, 5 (= F/cp 
where F is the local ionizing flux, c is the speed of light, 
and p is the gas pressure), and are valid in the optically 
thin, low-density regime. S is effectively a measure of the 
ratio of radiation pressure to the gas pressure, and is 
2.3P ra d/P gas for fully ionized gas of cosmic abundance. 
The heating rates include Compton and photoionization 
heating of all the ionization stages of hydrogen, helium, 
and some trace metals (particularly C, N, O, Si, S, and 
Fe). The cooling includes collisional excitation, recombina- 
tion cooling, Compton cooling, collisional ionization, and 
free-free l osses. The gas is assumed to be free o f dust (Laor 
& Draine 1993 ; but see review by Oster brock 1993] for an 
excellent summary of our current understanding of dust 
in AGN). The rates were calculated using the CLOUDY 
photoionization code, its standard AGN spectrum, and 
solar abundances. Whilst abundances in AGN remain a 
very contentious issue, there is considerable evidence that 
the nuclear regions are not metal poor, even at high red- 
shift, and also little evidence that the metallicity of the 
BELR changes with redshift (see Artymowicz et al. 1993). 
Though it seems certain that nitrogen is overabundant by 
factors of ~ 2 — 9, particularly for high-redshift sources 
(e.g. Hamann & Ferland 1992), most theoretical work has 
been based on the assumption that the gas is of solar 
composition: we also apply this assumption. More details 
of the heating and cooling rates can be found in Woods et 
al. ( |l99"6l ). 

In Fig. we show the thermal equilibrium curve 
for the assumed AGN spectrum. At low temperatures 
photoionization heating and cooling due to line excita- 
tion and recombination are in near balance. At large 
temperatures, the equilibrium arises from a balance of 
Compton heating and cooling. The solutions with p ositiv e 
slope, i.e. dT/dE > 0, are thermally stable (Field 1965 ). 
If the slope is negative, solutions can still be thermally 
stable if the gas cools isochorically, whilst they are ther- 
mally unstable if the gas cools isobarically. Note that 
there is a small range of ionization parameters for which 
there are stable solutions at intermediate temperatures 
(T ~ 10 6 K). The exact shape of this part of the ther- 
mal equilibrium curve is poorly known and most likely 
varies substantially from source to source, since it is a 
complicated function of the irradiating spectrum, the as- 
sumed abundances and thermal processes (cf. Krolik et 
al. |1981| ). For a given object it is entirely possible that 
there are no multi-valued equilibrium temperatures for 
any E. Therefore, in our discussion of the following results 
we do not assign too much importance to the exact be- 
haviour of the simulations in this part of parameter space 
since we are not modelling a specific object. In Sec. y 
we further show that the cooling gas often does not pass 
through this region of parameter space. 

To obtain cool gas in thermal equilibrium we require 
ionization parameters 3 < 10. As noted by Perry & Dyson 
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Fig. 2. Comparison of evolved radial profiles of density and temperature rescaled against the similarity solution (solid). 
The times are measured by the homologous expansion, with t = 1 corresponding to the initial mapping. The initial 
number of grid points in the region of shocked ejecta was four. The curves are not coincident because of start-up 
errors, finite numerical resolution, and unavoidable diffusion and dissipation in the numerical scheme. As the remnant 
expands the effective resolution of the calculation increases, and the profiles converge towards the original similarity 
solution. 



( |l985|) , shocked gas cooled back to equilibrium can have 
a value of 5 much lower than its pre-shock value. This 
is because 3 does not change if the gas cools isobarically 
and the post-shock pressure is much greater than the pre- 
shock value. Therefore strong shocks can create conditions 
for the gas to cool to temperatures much lower than the 
surrounding ambient temperature. 

The ionization parameter of the QSO ISM can be 
written as 



E w = 2 x 10 b 



Lboi n w r% c T 8 ' 



(7) 



where Lj „ and L\, \ are respectively the ionizing and 
bolomctric luminosity of the central source, L47 is the 
bolometric luminosity of the central source in units of 



10 ergs 



is the number density of the wind, 



is the distance of the shock from the central source in 



units of parsecs, and T 8 is the temperature of the wind in 
units of 10 8 K. 

The ionization parameter of shocked gas cooled back 
to equilibrium with the radiation field, S s fck, is given by 



1 — >shk 






(8) 



where u> is the velocity of the pre-shock gas relative to 
the shock in units of 0.01c (cf. Perry & Dyson 1985). For 
supernovae in the central regions of QSO, this velocity 
is a combination of the velocity of the blast wave, 
the velocity of the ambient medium, and the Keplerian 
velocity of the progenitor star. The QSO medium has 
characteristic velocities of v W j n d ~ 1, 000 — 3, OOOkms -1 
(u = 0.3 - 1.0, Williams et al. |l999|) , whilst the expan- 



sion velocities of young ejecta-dominated SNRs satisfying 
the Chevalier-Nadyozhin similarity solution can be greater 
than 10, OOOkms -1 . Therefore, it is not difficult to obtain 
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Fig. 3. Thermal equilibrium curve for the standard AGN 



spectrum in CLOUDY (see Woods et al. 1996). The sym- 
bols refer to the thermal evolution of the cool region 
in Models A (diamonds), B (triangles) and C (squares) 
(Fig. ||, Fig. U and Fig. ^ respectively). 



E shk < 1 for which T eq ~ 10 4 K (see Fig. g). The crucial 
question is whether the shocked gas remains at high pres- 
sures long enough to cool from its post-shock temperature 
to T ~ 10 4 K. This requires that the dynamical timescale 
for the SNR is longer than the cooling timescale of the 
post-shock gas. 



3. Results 

3.1. Adiabatic evolution 

We initially investigate the adiabatic evolution of a SNR 
surrounded by a static, constant density medium with 
n w = 10 6 cm" 3 and T w = 1.33 x 10 7 K. Whilst this 
value of n w is marginally optically thick to X-rays (at 
the Fe K-shell edge) and to Compton scattering for dis- 
tances greater than sw 0.5 pc, it is of no concern pro- 
vided that the density falls with radius on similar scales. 
However, as our results show, the radius of the supernova 
remnant remains significantly smaller than this scale. We 
note also that our chosen ambient density is similar to 
that found by Williams et al. (1999). In the central re- 
gions of AGN, the ISM density is primarily determined 
by mass-loss from the central stellar cluster through winds 
and supernova explosions. Collisions and tidal disruptions 
are negligible except in the central regions of the densest 
clusters. Compared to the central cluster, the mass flux 
from the external ISM (e.g. flow into the nucleus from a 
galactic bar) is unlikely to provide a significant fraction 
of the fuelling requirements of the AGN as a continuous 



source (Shlosman et al. 1990). Our value of T w is also 
characteristic of the Compton temperature in a hard QSO 
radiation field. 

The assumption of a static ambient medium keeps the 
problem in its simplest form. For the unshocked ejecta 
we set T = 10 4 K. The Chevalier-Nadyozhin similarity 
solution is derived on the assumption of negligible thermal 
energy in the unshocked ejecta and ambient medium, and 
our specified temperature implies thermal pressures orders 
of magnitude below the relevant ram pressure. Whilst this 
also implies exceedingly high absorption columns through 
the ejecta core, the covering fraction is very small. 

In Fig. [| the radius and speed of the contact 
discontinuity are shown as functions of time. As discussed 
earlier, the solution remains self-similar until the reverse 
shock propagates to the edge of the core: this occurs at 
(«9 yr, and is reflected in the structure of the curves in 
Fig. |. 

Although we have not included heating and cooling 
terms in this calculation, it is instructive to assume a 
radiation field and to plot the evolution of the ionization 
parameter of the equilibrium temperature of the shocked 
gas. This gives an impression of whether the shocked gas is 
likely to heat or cool if these terms are included. Assuming 
that the remnant is 0.33 pc distant from a central ioniz- 
ing source of L ion = 10 47 ergs _1 , we show in Fig. [| the 
time evolution of this parameter in the shocked ambient 
medium. This position was chosen because, with the as- 
sumption of constant ionizing flux, 5 tracks the inverse 
of p and the profile of p is flatter at this point. Fig. |^ 
shows that 5 initially evolves as i 1 / 2 , but increases more 
rapidly once the flow diverges from self-similarity. This is 
due to a large increase in the volume of the shocked re- 
gion (as the reverse shock begins to move back towards the 
centre of the remnant in the Lagrangian frame) at roughly 
constant total energy, which leads to a reduction in the in- 
ternal energy per unit volume and hence pressure, and the 
corresponding increase in S. It is therefore apparent that 
the formation of cool clouds in a SNR shock is favoured 
at early times, when the pressure of the shocked region is 
high and S is low. In particular, since S rapidly increases 
once the ejecta core reaches the reverse shock, a cool re- 
gion has the best opportunity to form before this point. 
Therefore we equate the dynamical timescale of the rem- 
nant, tdyni with the time at which the ejecta core reaches 
the reverse shock, and require that the shocked gas has 
a cooling time, t coo i, shorter than this. We derive expres- 
sions for tdyn and t CO oi in Appendix A. 

We also require that the ionization parameter 5 of the 
post-shock gas at t = td yn is low enough to allow the gas to 
cool to T w 2 x 10 4 K. Cold gas will not form if S > 2 cr it 
(where "E, C rit ^ 30 for the AGN spectrum used to generate 
the thermal equilibrium curve in Fig. y) whether or not 

tcool 



< 



tdyn- We derive an expression for this condition 
also in Appendix A. Hence for cold gas to form we need 
to satisfy t coo i < t dyn and E) tdyn < H C rit- Both of these 
conditions are evaluated for each of the following models. 
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Fig. 4. Evolution of the radius (solid) and speed (dots) of 
the contact discontinuity of an adiabatic SNR as a func- 
tion of time (arbitrary units). At t = 1 yr, their respective 
values are 1.0264 and 0.7698. The SNR evolves in a self- 
similar fashion with r ex £ 3 ' 4 and v oc £ -1 ' 4 until the 
reverse shock reaches the radius of the ejecta core. This 
occurs at t « 9 yr. Also shown is the time evolution of 5 
(dashes), which at t = 1 yr has a value of 1.7. Prior to 
t m 9 yr it evolves as t 1 ' 2 , but increases in value much 
more rapidly for t > 13 yr. 



3.2. Radiative evolution 

We now include the full heating and cooling rates in all 
of our remaining simulations. We start the evolution on 
the hydrodynamic grid at a time early enough for the 
dominant cooling of the shocked gas to be adiabatic ex- 
pansion. As the remnant expands, radiative cooling grad- 
ually increases in importance, and a smooth transition 
from the adiabatic self-similar solution into the radiative 
regime occurs. Because the post-shock temperatures can 
be very high early on (T > 10 9 K), we extrapolate the 
heating and cooling functions used by Woods et al. (1996) 
with a second order polynomial. This should be a rea- 
sonable approximation to the true rates because over the 
range 10 8 < T < 10 9 K, cooling and heating are domi- 
nated by Compton or bremmstrahlung processes which are 
smoothly varying functions of T. Furthermore, since the 
cooling of the remnant is initially dominated by adiabatic 
expansion, the inaccuracies in the correct rates from this 
extrapolation are small. Finally, as the remnant expands, 
the pre-shock ejecta are slowly heated from T — 10 4 K to 
~ 10 6 K. We assume again that the remnant is 0.33 pc dis- 



Table 1. Parameters for the models considered in this 
paper. All models have the same ionization parameter and 
temperature for the ambient medium (S « 150, T w — 
1.33 x 10 7 K). 



Model 



(cm- 3 ) (kms" 1 ) 



F 

- 1 ton 

(erg cm -2 s~ 



Cool 
Regions 



A 
B 

C 



Iff 
10 4 
10 4 





3000 



7.67 x 10 a 
7.67 x 10 7 
7.67 x 10 7 



Y 

N 
Y 



tant from a central source with L; 



10 47 ergs 4 . 



The 



parameters for this model (Model A) are listed in Table [y. 



In Fig. H the evolution of the shocked region is shown. 
We find that the shocked ejecta substantially cool with 
heating and cooling rates included in the calculation. By 
t = 1.0 yr a cool region has formed with a temperature 
less than that of the undisturbed ambient medium. At 
t = 1.4 yr it has a temperature of 7.0 x 10 5 K. Divergence 
from the adiabatic self-similar solution can be seen in 
the time evolution of the T, S, and p profiles. In the 
temperature profiles, this is first manifested as a change 
in the slope of the region of shocked ambient material 
from dT/dr > to dT/dr < 0. The formation of the 
cool region at later times is of course a much larger di- 
vergence from the self-similar solution. In the density pro- 
files, the formation of the cool region is apparent as a 
sharp coincident growth in density. The profiles of ioniza- 
tion parameter show the gradual temporal increase in 5 
expected from the evolution of the self-similar solution. 
However, at t = 1.4 yr, the value of S associated with the 
cool material is substantially larger than the value of the 
immediate surroundings. This is due to the radiative losses 
becoming so high that cooling no longer takes place isobar- 
ically. In the limit that the cooling timescale, t c , is much 
less than the appropriate dynamical timescale, td, the cold 
gas would cool isochorically. Here td is the timescale for 
the hot post-shock gas to respond to the rapid depressur- 
ization of this material, and is of the order of the timescale 
for collapse of the reverse shock, s=s l/(v c d — v rs ), where I 
is the length scale of the shocked ejecta, and v c d (v rs ) is 
the velocity of the contact discontinuity (reverse shock). 
The thermal parameters of the cool region at t = 0.5, 1.0, 
and 1.4 yr are plotted as crosses in Fig. ||. 

Beyond t = 1.4 yr, as the cloud cools further, its ther- 
mal energy as a fraction of its total energy (thermal plus 
kinetic) approaches the round-off error of our calculation 
(~ 10 -4 ), and we cannot follow its evolution past this 
time. There is also the well-known problem of the 'eat- 
ing away' of the edges of hot material. This affects all 
numerical schemes, and results from the numerics smear- 
ing the temperature gradient and placing cells at inter- 
mediate temperatures, which then undergo high cooling 
(often the line-cooling bump at T ~ 10 5 K). In this fash- 
ion cool regions can rapidly grow as hot regions are 'eaten 
away'. 

However, it is clear that the gas will continue to cool, 
and eventually reach T w 2 x 10 4 K, as the following 
argument demonstrates. With t c << td the gas cools iso- 
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chorically so p oc T oc 1/H. At t = 1.4 yr, the cloud is at a 
temperature T = 7.0 x 10 5 K and an ionization parameter 
H = 10.8. It therefore moves approximately along a line of 
slope -1 in Fig. 0, and encounters the thermal equilibrium 
curve at S « 40 and T w 2 x 10 5 K. This point is ther- 
mally stable and the rapid decrease in the temperature 
is brought to a halt. The temperature of the cool region, 
T c , then remains relatively constant, and ajusts only for 
changes in the value of the ionization parameter, S c , which 
occur on a sound-crossing timescale as the surrounding 
post-shock material gradually repressurizes the cool gas 
to its level. The sound crossing time within the shocked 
ejecta is short compared to the dynamical time of the rem- 
nant, so S c — > S ss , where S ss is the ionization parameter 
of the surrounding shock material (S ss w 2 — 3 in the 
shocked ambient gas). There is therefore no immediate 
danger of S c increasing to the point that the cooled gas 
begins to reheat towards the Compton temperature. Thus 
S c approaches S ss and the cooled region moves along the 
thermal equilibrium curve to the left. The density of the 
region, n c , increases during this process until 



T r 



(9) 



where n c i (T c > ) is the density (temperature) of the cooled 
material at the end of our simulation. From the results in 
Fig. |[ we obtain n c ~ few xlO 10 cm -3 , in good agreement 
with inferred values of the density of the HIL region from 
observations. We can also estimate the resultant hydrogen 
column of the region: at t = 1.4 yr, its thickness is approx- 
imately 2 x 10 13 cm (the numerical resolution of our calcu- 
lation is about half this), and its density is rs 3 x 10 8 cm 
which for solar abundances results in a column density of 
w 10 22 cm~ 2 , again in agreement with observations of 
the HILs. The radius of the remnant is « 7.5 x 10~ 3 pc 
at this stage, whilst the expansion speed of the contact 
discontinuity (and thus the velocity of the cool region) is 
w 4 x 10 8 cms -1 . 

Observations of the HILs have revealed that the op- 
tical depth is less than but of order unity, and that 
Nh ~ 10 22 cm -2 . The observed range of the ionization 
paramet er of the HIL region is 0.3 < S < 2 (Kwan & 
Krolick 1981 ), which is somewhat higher than that in- 
ferred for the LIL region. The systematic blue-shift of 
the HILs with respect to the LILs is interpreted as the 
HILs forming in a region with bulk outflow (possibly in- 
flow) and that the red-shifted emission is obscured. Line 
widths, which reflect local gas velocities, are typically sev- 
eral 10 3 kms _1 . Hence our model results are in harmony 
with the temperature, ionization parameter, column den- 
sity, and velocity dispersion of the observed HIL BELR 
clouds. 

We now check to see if our results are in agree- 
ment with the two conditions which must be satisfied 
for the shocked gas to cool down to T rs 2 x 10 4 K, 
namely t cool < t dyn and S) tdyn < E crit . For n = 12, 
S = 0, E = 10 51 ergs, M ej = 10 M®, R 2 /R c = 0.974, 



tdyn = 7 yr. Assuming solar abundance for the ejecta, 
(i s = 0.61, and t coo i = 3 yr. Hence t coo i < t dyn as re- 
quired. The fact that the cooled region first forms at 
t » 1.4 yr indicates that we are underestimating the cool- 
ing rate at lower temperatures where line-cooling is dom- 
inant. For these parameters, we further obtain from Eq. U\ 

£shk)t dvn = 6.3, 



V r 



3610 kms , and from Eq. A. 13 



dyn 



and thus Model A also satisfies the requirement that 
^■shk)t d „ ^ 20. In reality the ejecta will be enriched in 
metals, which will reduce the above estimate of t coo i. 



3.3. Exploration of associated parameter space 



3.3.1. Variation of n„ 



r pc , and L z 



Since n w , r pc , and Li on all influence S s ^fc it is impossible 
to discuss the effect of a variation in one without also con- 
sidering the impact of the others. This inevitably leads to 
some complexities but we attempt to keep our discussion 
as simple as possible. 

We can first consider the effect of varying the assumed 
ambient density. Since the equilibrium temperature is a 
function of 5, we wish to keep this value constant for 
our initial investigation. As E w oc Li on /(n w ri c ), this im- 
plies a variation in flux from the central engine, which 
can obviously be interpreted as a variation in the ioniz- 
ing luminosity L.; on with r pc fixed, or as a variation in the 
distance r pc between the remnant and the central engine 
with Li on fixed, or a suitable variation in both Li on and 



3 CT 



A = 0.19 and p w = 10 8 g cm , we obtain from Eq. A. 5 



In Fig. g we show the results for a simulation with 
~ 150 (matching the earlier simulation of Model A 
shown in Fig. H) but with a lower ambient density 
[n w = 10 4 cm" 3 ) and AGN flux (e.g. L ion = lO^ergs" 1 
with the remnant 1 pc distant). We call this Model B. As 
a result of the lower ambient density, it has a longer evolu- 
tion time (tdyn = 32 yrs) and the remnant expands further 
than in Model A. Relative to Model A, the lower flux from 
the central engine also increases the cooling timescale of 
the shocked gas (t coo i = 300 yrs). Although the ionzia- 
tion parameter of the shocked gas at t = tdyn is less than 
S c „ t (S) tdH „ = 6.4, as for Model A), because t coo i > t dyn 
we do not expect any of the shocked material to cool to 
T w 2 x 10 4 K, 

Once the reverse shock reaches the core radius at 
t = 32 yr, the solution swiftly diverges from a self-similar 
form (despite the cooling of the post-shock material, many 
aspects such as the radius and velocity of the shocks are 
not too far removed from such an evolutionary form prior 
to this point). At this time the coolest shocked gas is at 
T w 7 x 10 6 K. The subsequent acceleration of the re- 
verse shock towards the centre of the remnant rapidly 
de-pressurizes the shocked region, increasing the ioniza- 
tion parameter of the shocked gas, and ultimately leading 
to its reheating towards the Compton temperature. Thus, 
by simply adjusting the density of the ambient medium, 
whilst keeping all other parameters constant, we have 
shown that a cool region may not form if the effective 
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Fig. 5. The evolution of the remnant of a type-II supernova for a constant density medium with n w = 10 6 cm -3 and 
'Bw « 150 (Model A). The medium is stationary with respect to the remnant, and the entire area is immersed in 
a radiation field appropriate to an AGN of ionizing luminosity 10 47 ergs _1 with the central engine at a distance of 
0.33 pc. The panels show from top to bottom the evolution of the temperature, ionization parameter, and density of 
the region of hot shocked gas. This is bounded on its right edge by the forward shock propagating into the ambient 
medium and on its left edge by the reverse shock. From left to right the corresponding times are 0.2 yr, 0.5 yr, 1.0 yr 
and 1.4 yr. The formation of a cooled region of gas can clearly be seen. 



ram pressure of the ambient medium is not sufficiently 
high. 

We consider whether there is also a mechanism which 
prevents the formation of cool regions in SNRs once the 
density of the ambient material increases passes some 
limit. This question can be addressed in the following way. 
As n w increases, Li on /r"i c must similarly increase to main- 
tain a given value for B w . For a given Li on this means a re- 
duction in Tp C (which might be consistent with an increase 
in n w ). A higher value of n w leads to a faster evolution 
of the remnant, which given that the cooling timescale 
of the shocked gas also shortens, does not prevent the 
possibility of cool regions forming. However, physically it 
does mean that the lifetime of the cool regions is also very 
short, since the evolved time between their formation and 
their destruction, which presumably occurs shortly after 
the interaction of the ejecta core with the reverse shock, is 
short too. Hence at some value of n w , the contribution to 
the BELR emission will be on such short timescalcs that 



it will become an unimportant component of the overall 
emission. 

We note that it is unlikely that the ejecta density actu- 
ally has such a sharp cut in its gradient. However, once the 
ejecta passing through the reverse shock starts deviating 
from an r~ 12 profile the volume of the shocked region will 
start expanding with an inevitable rise in the ionization 
parameter. This would occur at an earlier time relative to 
the sharp cut-off case and would tend to reduce the con- 
tribution of mass to the BELR by this model. Assuming 
a constant value of B w and Li on , this argument also in- 
troduces a lower bound to the value of r pc at which a cool 
region can form. 

Conversely, for a given r pci an increase in n w and con- 
stant Sk, means a corresponding increase in Li on . In a sim- 
ilar way to the above argument, this introduces an upper 
bound to the value of Li on at which a cool region can 
form. These conditions together place upper bounds on 
the values of n r and v r inferred from observations. 
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Fig. 6. As Fig. a (Model A) but for a constant density medium with n w = 10 4 cm~ 3 , an ionizing AGN luminosity 
of 10 46 ergs _1 , and a distance of 1 pc to the central engine (Model B). E w is again w 150. From left to right the 
corresponding times are 20 yr, 50 yr, 120 yr and 180 yr. In contrast to Model A, cool clouds have no opportunity to 
form since the interaction of the ejecta core with the reverse shock leads to a rapid depressurization of the shocked 
gas. 



3.3.2. Variation of E w 

We now consider the effect of a different wind ionization 
parameter on the possible formation of a cool BELR. In 
our simulations so far we have deliberately specified a 
value of E w which is high enough for the ambient medium 
to exist at the Compton temperature, but which is at the 
same time almost as low as we could make it, providing us 
with the best opportunity to form a cool region (through 
a low value of E s hk cf. Eq. ||). If 5^ was increased in value 
(and given that T w is almost constant for large values of 
E w ) we would require a higher relative velocity, u>, to ob- 
tain a given value of E s hk (cf. Eq. 0). For a given remnant 
age and n w , and assuming a static medium, a higher value 
of E w gives a higher value of E s hk, and since E s hk grad- 
ually increases with time as the expansion velocity of the 
remnant slows, this effectively reduces the 'window' of op- 
portunity where 'E.shk is low enough for T eq ~ 10 4 K. This 
in turn places a limit on the maximum cooling timescale 
for the shocked gas, which for a given value of E w , feeds 
back into a lower limit on the ambient density, n w ^ m { n . If 



n w < n WyTn i ni the shocked gas will not have enough time 
to cool to ~ 10 4 K before E s hk becomes too large. 

3.3.3. Variation of ui 

If the medium surrounding the remant has some bulk 
velocity of its own, it is possible to obtain the same value 
of E s hk at the same remnant age, for a higher value of E w . 
On the side of the remnant facing the oncoming 'AGN 
wind' the expansion of the remnant is inhibited, leading 
to an increased value for v re i (cf. Eqs. |A.8| and A.G) or 
alternatively w (cf. Eq. @). In turn this reduces the critical 
density of the ambient medium needed for a cool region to 
form, n w ^ min . In Fig. |we show the results for Model C, a 
remnant of a type-II supernova expanding into a constant 
density medium with n w — 10 4 cm" 3 , E w ss 150 , and 
a flow velocity of 3, 000 km s" 1 (cf. Williams et al. |l999| ) 
towards the remnant. These results are indicative of the 
remnant structure on the upstream side (since the shocked 
region is thin with respect to the radius of the remnant). 
The entire area is again immersed in a radiation field ap- 
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propriate to a QSO of luminosity 10 46 ergs~ 4 with the 
central engine at a distance of 1 pc. Fig. R shows that 
this time a cool region can form. In this case the ioniza- 
tion parameter of the cool gas is expected to be ~ 0.5 
and its density to be n « 10 9 cm -3 . Its column density is 
w 6 x 10 21 cm~ 2 and its speed is « 2580kms _1 . These 
values are again in good agreement with observations. 



3.3.4. Constraints on the parameter range 

The previous sections show that there are two broad con- 
straints on the parameter range required for cool regions 
to form in SNRs. First, the post-shock ionization param- 
eter needs to be low enough for the equilibrium tempera- 
ture of this region to approach w 10 4 K. This requires that 
the SNR evolves in a dense medium (Model A (Fig. H) vs. 
Model B (Fig. o)) but at the same time is not too close 
to the central engine (cf. Eq. A.13J ). Whilst tending to de- 
crease the evolutionary timescale of the remnant, a higher 
effective confining pressure increases the density of the 
shocked gas leading to shorter cooling timescales. Thus 
although a high ambient density does not inhibit the for- 
mation of a cool region, at some level it results in such 
rapid evolution that these regions are destroyed without 
making any substantial contribution to the HIL emission. 

The distance to the central engin e inf luences the value 
of S s hfe through the ionzing flux (Eq. |A.8| ) Therefore there 
is an upper limit to 3^ after which not even SNe in the 
densest environments will be able to produce cool regions. 
If the ambient medium is in motion relative to the rem- 
nant, the parameter space over which cool regions can 
form is extended (Model B (Fig. |) vs. Model C (Fig. @)). 

The second constraint is that the ambient medium 
needs to be dense enough for the cooling timescale of the 
shocked gas to be less than the age of the remnant at 
the time when the reverse shock reaches ejecta which do 
not have a strongly radially dependent density. If this is 
not the case no cool regions can form, as illustrated by 
Model B (Fig. |). 



4. Discussion 

The results in Sec. || demonstrate that it is possible to 
cool shocked supernova ejecta down to T ~ 10 4 K in 
the inner regions of a QSO. Although our results differ 
from the original proposals of Dyson & Perry ( 1982) and 
Perry & Dyson ( |1985| ), which were for the shocked am- 
bient medium to cool, the resulting cool gas nevertheless 
has properties (densities, column depths, velocities and 
ionization parameters) compatible with those inferred for 
gas emitting the high ionization lines in QSOs. We have 
additionally shown from three separate calculations that 
this model is able to cover the parameter range associ- 
ated with observations of BELRs, and have highlighted 
scenarios which are outside of this range and preclude the 
formation of cool regions. 

We now compare our results to some of the obser- 
vational correlations that have been determined. Korista 



et al. (1995) noted that the blue and red wings respond 
faster to continuum changes than the core. If one as- 
sumes that the expansion velocity of the supernova rem- 
nant is the dominant velocity component, this observation 
can be explained with our model in the following way. 
Remnants which are surrounded by the highest ambient 
densities evolve most rapidly, and therefore form cool re- 
gions with the highest expansion velocities. If the ambient 
density increases towards the central engine, there will be 
a tendency for remnants in these regions to form cool re- 
gions with, on average, higher velocities. Being close to the 
central engine, these regions respond more rapidly to any 
continuum changes, as the observed correlation requires. 
Even if v exp is not the dominant component, the above 
argument implies that our assumed model nevertheless 
provides a statistical bias in this direction. 

There now appears to be conclusive evidence that the 
higher ionization lines respond faster than lower ioniza- 
tion lines to changes in the continuum variability (see ref- 
erences in Fromerth & Melia 2001). If the LIL are formed 
in an accretion disc close to the central engine, this implies 
that the thermal timescale in the disc is longer than the 
corresponding timescale for the HIL. Alternatively the de- 
lay could be explained if the bulk of the LIL heating is due 
to back-scattered X-rays from the gas responsible for the 
HILs and the general ambient medium (see Collin-Souffrin 
et al. 1988). Therefore, our model is consistent with these 
observations. 

Our calculations show that it is the shocked ejecta 
which cool to form the HILs, and that once the ejecta core 
reaches the reverse shock, the shocked region rapidly de- 
pressurizes. Since this presumably leads to the destruction 
of any cool regions, it seems sensible to suppose that for 
the remnant parameters we have chosen, an upper limit 
to the HIL mass would be 2.5 M Q per remnant. With 
further appropriate assumptions we can then estimate the 
required supernova rate needed to sustain the mass in the 
gas responsible for the HIL. With the assumption that 
each remnant can eventually cool ~ 2.5 M of gas, and 
that the lifetime of the cool gas in the remnant is of order 
10 yr, to sustain ~ 25 M of gas in the HIL would require 
a supernova rate of ~ 2 yr -1 . Although our estimates are 
very uncertain, this rate is acceptable: certainly for high 
luminosity QSOs, supernova rates ~ 10 yr -1 are conceiv- 



able (Terlevich et al. 1992J ). A less steep density distribu- 
tion of the ejecta envelope (if n < 12 more mass would be 
available per SN) , and the suppression of dust formation in 
intermediate-mass AGB stars in the BELR region (which 
may reduce the minimum zero-age main sequence mass re- 
quired for supernovae; Hartquist et al. 1998D are two addi- 
tional possibilities which could further reduce our estimate 
of the required supernova rate. Additionally, it is possible 
that conditions exist for the swept-up ambient material to 
also cool to temperatures appro pria te for the HIL, as the 
cooling time estimated by Eq. A. 7 is temperature inde- 
pendent. The fact that the reverse shocked material cools 
first in our simulations is simply due to enhanced cool- 
ing by line emission starting earlier due to its lower initial 
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Fig. 7. As Fig. o (Model B) but in this simulation (Model C) the ambient medium has a flow velocity towards the 
remnant of 3, OOOkms -1 . From left to right the corresponding times are 5 yr, 10 yr, 12 yr and 12.5 yr. In contrast 
to the results for a stationary medium, cool clouds can form in this case. Note also the more rapid evolution of the 
remnant with respect to the case shown in Fig. o. 



temperature. In Model A where t, 



cool 



<t, 



dyn 



we therefore 



might expect the swept-up ambient material to form HIL 
gas as well. At t = td yn , the forward shock is estimated 
to be at a radius of 0.03 pc and to have swept up 1.6 M 
of ambient gas. This gas could therefore be a significant 
contributer to the total HIL emission. 

In this first paper we have been solely interested in 
the question of whether cool regions could form from su- 
pernova shocks given that they are bathed in the hard 
radiation flux of the central engine. For simplicity we re- 
stricted our modelling to the simplest ID approach, and 
assumed solar metallicites in our calculations. The fact 
that the ejecta (and also the swept up medium) may be 
responsible for the HIL emission warrants a careful consid- 
eration in future models. We will also perform calculations 
on 2D axisymmetric hydrodynamic grids. 

Another future goal is the inclusion of the effect of the 
QSO radiation field on the dynamics of the remnant. At 
high temperatures, the effective cross-section of gas is the 
Thomson scattering cross-section, ut- Once the gas begins 
to cool, the effective cross-section, a, increases, enh ancin g 
the radiative driving. As noted by Williams et al. ( 1999J) , 



a/ar increases to roughly 2 at 10 6 K, 40 at 10 5 K, and 
to > 10 4 when the gas has fully cooled (see also Arav 



& Li ( |1994| ) and Arav et al. ( |1994| )). For gas cooler than 
approximately 5 x 10 4 K, it is reasonable to take the ra- 



diati ve acc eleration as g — 7 x 10 14 p (Roser 1979] ; Dyson 



et al. 1981 ). We expect that substantial radiative driving 
will critically alter the dynamics of the interaction be- 
tween the supernova ejecta and the ambient medium (e.g. 
Falle et al.|l98l| Williams |2000|). 



Although we have chosen to model the interaction of 
supernova ejecta with the ambient medium, it is possible 
that a wind from a group of early-type stars may also pro- 
vide the necessary conditions for the formation of cool re- 
gions. This interaction may be more relevant in the nuclei 
of Seyfert galaxies, since supernova explosions will evac- 
uate all but the most tightly bound gas in them (Perry 
& Dyson 1985). Finally it is clear from our models that 
whilst the supernova-QSO wind interaction is conceptu- 
ally simple, the BELR is likely to be a very complicated 
region in practice. 
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Appendix A: Derivation of td yn and H 



slik 



From Chevalier (1982), the radius of the contact 



discontinuity for a supernova remnant expanding 
into a stationary, constant density ambient medium 
(p w = liuj^niH, where jj, w is the average mass per 
particle), is 



Re 



Ag r 



n w (J, w mH 



1/,. 



(n— 3)/n 



(A.l) 



where A is a constant determined by the steepness of the 
ejecta profile (n). The radius of the reverse shock is 



Ri — ~~^~Rc 
R c 



(A.2) 



where R2/R c is also a function of n. For n — 12, 
Chevalier's Table 1 gives A = 0.19 and R 2 /R c = 0.974. 
The velocity of the pre-shock ejecta at the reverse shock 

is 



v ej ) R2 = R 2 /t. 



(A.3) 



The dynamical time (i.e. the time at which the ejecta core 
reaches the position of the reverse shock) can be obtained 
from the ratio 



I'r 



tdyn 



-n/3 



(A.4) 
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where v e j)n 2 is the velocity of the pre-shock ejecta at the Combining Eqs. [A.8| , A. 9 and A. 12, and noting that as 
reverse shock at time t. Substituting Eqs. nL H and A. 3 we it is the reverse shock which cools, n s = An w p2Jp\ = 
obtain 



tdyr. 



A 



1/3 



Ih 



^/3 



47t(n - 5)n w p w nin 



[(3-S)(n-3)M ej ]-s/v 



(A.5) 



Since the shocked gas is initially at very high temper- 
atures, we make the assumption that the dominant cool- 
ing during this interval is inverse Compton. The radiative 
losses (erg cm" 3 s _1 ) are given by 



E 



2.2 x 10 7 ^T 8 p 



(A.6) 



4pwP2/(piHw hih), we obtain 
1 p\ Li, 



^shk )td 



J ionA7 



27 p2 r% c n w p w m H I [l - ( Ii ^)] i 
Values for P1/P2 are given in Chevalier (|1982) 



.(A.13) 



(Perry & Dyson 1985). The cooling timescale is indepen- 
dent of density and temperature, and is given by 



tcool = — ~ 18 

E PsLtf 



yrs, 



(A.7) 



where p s is the average mass per particle for the shocked 
material. For the post-shock gas to cool down to approxi- 
mately T = 2 x 10 4 K, we require t coo i < ti yn . 

We now derive an expression for the evolution of the 
ionization parameter of post-shock gas cooled back to 
equilibrium with the radiation field, E a . From Perry & 



Dyson (1985) we have 

— l-ti.nn 



^shk 



A-Kr 2 n s kT S ' 



2.0 x 10 



14 '-'ion, 47 



r pc n s l s 



(A. 



where n s is the number density of the shocked material, T s 
is the immediate post-shock temperature of the material 
and is given by 



3 /U s m H 2 

Ts= 16^^' 



(A.9) 



and v r is the velocity of the pre-shock ejecta in the frame 
of the reverse shock. The velocity of the reverse shock is 



VR 2 



R 2 
R, 



■Vcd 



3i? 2 



R c 



Ag n 



n w fi w rim 



l/r, 



-3/n 



(A.10) 



so using Eq. A. 3 we obtain 



1- 



n-3\j R 2 ( A 

n )\ R c \4irn w iJ, w mTi(n — 6) 
[2(5 -S)(n-5)E}^)/^ 3/n 

[(3-<5)(n-3)M ej ](™- 5 )/ 2 » 

which at t = tdyn simplifies to 
' n — 3 



l/n 



v r )t 



1 



(A.ll) 



(A.12) 



